This notebook demonstrates how to perform phase and electrochemical assessments starting from a VASP calculation using Python Materials Genomics (pymatgen) and the Materials Project database (via the Materials API). These notebooks are described in detail in
Deng, Z.; Zhu, Z.; Chu, I.-H.; Ong, S. P. Data-Driven First-Principles Methods for the Study and Design of
Alkali Superionic Conductors. Chem. Mater. 2017, 29 (1), 281–288 DOI: 10.1021/acs.chemmater.6b02648.
If you find these notebooks useful and use the functionality demonstrated, please consider citing the above work.
Let's start by importing some modules and classes that we will be using.
In [1]:
%matplotlib inline
from pymatgen import MPRester, Composition, Element
from pymatgen.io.vasp import Vasprun
from pymatgen.phasediagram.maker import PhaseDiagram, CompoundPhaseDiagram
from pymatgen.phasediagram.analyzer import PDAnalyzer
from pymatgen.phasediagram.plotter import PDPlotter
from pymatgen.entries.computed_entries import ComputedEntry
from pymatgen.entries.compatibility import MaterialsProjectCompatibility
from pymatgen.util.plotting_utils import get_publication_quality_plot
import json
import re
import palettable
import matplotlib as mpl
In [2]:
vasprun = Vasprun("aimd_data/vasprun.xml.relax2.gz")
# include structure so proper correction can be applied for oxides and sulfides
entry = vasprun.get_computed_entry(inc_structure=True)
To construct the phase diagram, we need all entries in the Li-P-S-Cl chemical space. We will use the MPRester class to obtain these entries from the Materials Project via the Materials API.
In [3]:
rester = MPRester()
mp_entries = rester.get_entries_in_chemsys(["Li", "P", "S", "Cl"])
In addition to all the MP entries, here we also load the computed entries of O/S substituted Li-P-O tenary compounds.
In [4]:
with open("aimd_data/lpo_entries.json") as f:
lpo_data = json.load(f)
lpo_entries = [ComputedEntry.from_dict(d) for d in lpo_data]
Next, we need to combine all the entries and postprocess them using MaterialsProjectCompatibility. This postprocessing step corrects the energies to account for well-known DFT errors, e.g., in the sulfur binding energy.
In [5]:
compatibility = MaterialsProjectCompatibility()
entry = compatibility.process_entry(entry)
entries = compatibility.process_entries([entry] + mp_entries + lpo_entries)
In [6]:
pd = PhaseDiagram(entries)
plotter = PDPlotter(pd)
plotter.show()
We may observe from the above phase diagram that Li6PS5Cl is not a stable phase (red nodes) in the calculated 0K phase diagram.
The pseudo-ternary Li2S-P2S5-LiCl is constructed using the CompoundPhaseDiagram class.
In [7]:
cpd = CompoundPhaseDiagram(entries,
[Composition("P2S5"), Composition("Li2S"), Composition("LiCl")])
cplotter = PDPlotter(cpd, show_unstable=True)
cplotter.show()
In [8]:
analyzer = PDAnalyzer(pd)
ehull = analyzer.get_e_above_hull(entry)
print("The energy above hull of Li6PS5Cl is %.3f eV/atom." % ehull)
In [9]:
li_entries = [e for e in entries if e.composition.reduced_formula == "Li"]
uli0 = min(li_entries, key=lambda e: e.energy_per_atom).energy_per_atom
The PDAnalyzer class provides a quick way to plot the phase diagram at a particular composition (e.g., Li6PS5Cl) as a function of lithium chemical potential called get_element_profile.
In [10]:
el_profile = analyzer.get_element_profile(Element("Li"), entry.composition)
for i, d in enumerate(el_profile):
voltage = -(d["chempot"] - uli0)
print("Voltage: %s V" % voltage)
print(d["reaction"])
print("")
This element profile can be plotted as a Li evolution versus voltage using matplotlib as follows.
In [11]:
# Some matplotlib settings to improve the look of the plot.
mpl.rcParams['axes.linewidth']=3
mpl.rcParams['lines.markeredgewidth']=4
mpl.rcParams['lines.linewidth']=3
mpl.rcParams['lines.markersize']=15
mpl.rcParams['xtick.major.width']=3
mpl.rcParams['xtick.major.size']=8
mpl.rcParams['xtick.minor.width']=3
mpl.rcParams['xtick.minor.size']=4
mpl.rcParams['ytick.major.width']=3
mpl.rcParams['ytick.major.size']=8
mpl.rcParams['ytick.minor.width']=3
mpl.rcParams['ytick.minor.size']=4
# Plot of Li uptake per formula unit (f.u.) of Li6PS5Cl against voltage vs Li/Li+.
colors = palettable.colorbrewer.qualitative.Set1_9.mpl_colors
plt = get_publication_quality_plot(12, 8)
for i, d in enumerate(el_profile):
v = - (d["chempot"] - uli0)
if i != 0:
plt.plot([x2, x2], [y1, d["evolution"] / 4.0], 'k', linewidth=3)
x1 = v
y1 = d["evolution"] / 4.0
if i != len(el_profile) - 1:
x2 = - (el_profile[i + 1]["chempot"] - uli0)
else:
x2 = 5.0
if i in [0, 4, 5, 7]:
products = [re.sub(r"(\d+)", r"$_{\1}$", p.reduced_formula)
for p in d["reaction"].products if p.reduced_formula != "Li"]
plt.annotate(", ".join(products), xy=(v + 0.05, y1 + 0.3),
fontsize=24, color=colors[0])
plt.plot([x1, x2], [y1, y1], color=colors[0], linewidth=5)
else:
plt.plot([x1, x2], [y1, y1], 'k', linewidth=3)
plt.xlim((0, 4.0))
plt.ylim((-6, 10))
plt.xlabel("Voltage vs Li/Li$^+$ (V)")
plt.ylabel("Li uptake per f.u.")
plt.tight_layout()
In [ ]: